home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Cream of the Crop 22
/
Cream of the Crop 22.iso
/
math
/
ast53src.zip
/
PLACALC2.C
< prev
next >
Wrap
C/C++ Source or Header
|
1996-09-29
|
36KB
|
1,002 lines
/*
** Astrolog (Version 5.30) File: placalc2.c
**
** IMPORTANT NOTICE: The graphics database and chart display routines
** used in this program are Copyright (C) 1991-1996 by Walter D. Pullen
** (Astara@msn.com, http://www.magitech.com/~cruiser1/astrolog.htm).
** Permission is granted to freely use and distribute these routines
** provided one doesn't sell, restrict, or profit from them in any way.
** Modification is allowed provided these notices remain with any
** altered or edited versions of the program.
**
** The main planetary calculation routines used in this program have
** been Copyrighted and the core of this program is basically a
** conversion to C of the routines created by James Neely as listed in
** Michael Erlewine's 'Manual of Computer Programming for Astrologers',
** available from Matrix Software. The copyright gives us permission to
** use the routines for personal use but not to sell them or profit from
** them in any way.
**
** The PostScript code within the core graphics routines are programmed
** and Copyright (C) 1992-1993 by Brian D. Willoughby
** (brianw@sounds.wa.com). Conditions are identical to those above.
**
** The extended accurate ephemeris databases and formulas are from the
** calculation routines in the program "Placalc" and are programmed and
** Copyright (C) 1989,1991,1993 by Astrodienst AG and Alois Treindl
** (alois@azur.ch). The use of that source code is subject to
** regulations made by Astrodienst Zurich, and the code is not in the
** public domain. This copyright notice must not be changed or removed
** by any user of this program.
**
** Initial programming 8/28,30, 9/10,13,16,20,23, 10/3,6,7, 11/7,10,21/1991.
** X Window graphics initially programmed 10/23-29/1991.
** PostScript graphics initially programmed 11/29-30/1992.
** Last code change made 9/22/1996.
*/
#include "placalc.h"
#ifdef PLACALC
/*
** ---------------------------------------------------------------
** | Copyright Astrodienst AG and Alois Treindl, 1989,1991,1993 |
** | The use of this source code is subject to regulations made |
** | by Astrodienst Zurich. The code is NOT in the public domain.|
** | |
** | This copyright notice must not be changed or removed |
** | by any user of this program. |
** ---------------------------------------------------------------
**
** Important changes:
** 11-jun-93 revision 1.12: fixed error which affected Mercury between -2100
** and -3100 (it jumped wildly).
*/
#ifdef ASTROLOG
/* Given an object index and a Julian Day time, get its zodiac and */
/* declination position (planetary longitude and latitude) of the object */
/* and its velocity and distance from the Earth or Sun. This basically */
/* just calls the Placalc calculation function to actually do it, but as */
/* this is the one routine called from Astrolog, this is the one routine */
/* which has knowledge of and uses both Astrolog and Placalc definitions, */
/* and does things such as translation to Placalc indices and formats. */
bool FPlacalcPlanet(ind, jd, helio, obj, objalt, dir, space)
int ind, helio;
double jd;
real *obj, *objalt, *dir, *space;
{
int iobj, flag;
REAL8 jd_ad, rlng, rrad, rlat, rspeed;
if (ind <= oPlu) /* Convert Astrolog object index to Placalc index. */
iobj = ind-1;
else if (ind == oChi)
iobj = CHIRON;
else if (FBetween(ind, oCer, oVes))
iobj = ind - oCer + CERES;
else if (ind == oNod)
iobj = us.fTrueNode ? TRUE_NODE : MEAN_NODE;
else if (ind == oLil)
iobj = LILITH;
else
return fFalse;
jd_ad = jd - JUL_OFFSET;
flag = helio ? CALC_BIT_SPEED | CALC_BIT_HELIO : CALC_BIT_SPEED;
jd_ad += deltat(jd_ad);
if (calc(iobj, jd_ad, flag, &rlng, &rrad, &rlat, &rspeed) == OK) {
*obj = rlng;
*objalt = rlat;
*dir = rspeed;
*space = rrad;
return fTrue;
}
return fFalse;
}
#endif /* ASTROLOG */
/***********************************************************
** $Header$
**
** definition module for planetary elements
** and disturbation coefficients
** version HP-UX C for new version with stored outer planets
** 31-jul-88
** by Alois Treindl
**
** ---------------------------------------------------------------
** | Copyright Astrodienst Zurich AG and Alois Treindl, 1989. |
** | The use of this source code is subject to regulations made |
** | by Astrodienst Zurich. The code is NOT in the public domain.|
** | |
** | This copyright notice must not be changed or removed |
** | by any user of this program. |
** ---------------------------------------------------------------
**
***********************************************************/
/************************************************************
externally accessible globals, defined as extern in placalc.h
************************************************************/
REAL8 meanekl, ekl, nut;
struct elements el[MARS + 1];
/*
** In the elements degrees were kept as the units for the constants. This
** requires conversion to radians, when the actual calculations are performed.
** This approach is not the most efficient, but safer for development.
** Constant conversion could be done by writing all degree constants with
** value * DEGTORAD
*/
#define TIDAL_26 TRUE /* decide wheter to use new or old lunar tidal
term; a consistent system of delta t must be
used */
#define MOON_TEST_CORR FALSE /* to include more lunar terms in longitude */
REAL8 ekld[4] = {23.452294, -46.845, -.0059, 0.00181};
/* ecliptic with epoch1900, Ekd(0..3) in basic */
struct eledata pd[MARS + 1] = {
{/*earth*/ 1.00000023, 365.25636042, EPOCH1900,
99.696678, .9856473354, 1.089, 0,
101.220833, 6189.03, 1.63, 0.012,
0.01675104, -0.00004180, -0.000000126,
0, 0, 0, 0,
0, 0, 0},
/*
** note 29 June 88 by Alois: G.M.Clemence, Astronomical Journal
** vol.53,p. 178 (1948) gives a correction to the perihel motion
** of -4.78" T, giving 6184.25 for the linear Term above. We have
** not yet applied this correction. It has been used in APAE 22,4
** on the motion of mars and does make an official impression.
*/
{/*moon*/ 0.0025955307, 27.321661, EPOCH1900,
# if ! TIDAL_26
/*
** values from Improved Lunar Ephemeris, corresponding to tidal
** term -22.44"/cy and consistent with delta t ~ 29.949 T*T
*/
270.4341638, 13.176396526808121, -4.08, 0.0068,
# endif
# if TIDAL_26
/*
** new values from Morrison 1979, with tidal term -26"/cy as
** stated in A.E. 1986 onwards, consistent with delta t ~ 44.3 T*T
** correction: -1.54" + 2.33" T - 1.78" T*T
*/
270.4337361, 13.176396544528099, -5.86, 0.0068,
# endif
334.329556, 14648522.52, -37.17, -0.045,
0.054900489, 0, 0,
259.183275, -6962911.23, 7.48 , 0.008,
5.145388889, 0, 0},
{/*mercury*/ .3870986, 87.969252, EPOCH1900,
178.179078, 4.0923770233, 1.084, 0,
75.89969722, 5599.76, 1.061, 0,
0.20561421, .00002046, -.000000030,
47.145944444, 4266.75, .626, 0,
7.0028805555, 6.699, -.066},
{/*venus*/ .72333162, 224.700726, EPOCH1900,
342.767053, 1.6021687039, 1.1148, 0,
130.16383333, 5068.93, -3.515, 0,
0.00682069, -.00004774, .000000091,
75.7796472223,3239.46, 1.476, 0,
3.3936305555, 3.621, .0035},
{/*mars*/ 1.5236914620, 686.9296097, EPOCH1900,
/* These are the corrected elements by Ross */
293.74762778, .524071163814, 1.1184, 0,
334.21820278, 6626.73, .4675, -0.0043,
0.09331290, .000092064, -.000000077,
48.786441667, 2775.57, -.005, -0.0192,
1.85033333, -2.430, .0454}
};
struct sdat _sd [SDNUM] = {
114.50, 585.17493,
109.856, 191.39977,
148.031, 30.34583,
284.716, 12.21794,
114.508, 585.17656,
-0.56, 359.99213,
148.03, 30.34743,
284.72, 12.2196,
248.07, 1494.726615,
359.44, 359.993595,
109.86, 191.402867,
148.02, 30.348930,
114.503, 585.173715,
359.444, 359.989285,
148.021, 30.344620,
284.716, 12.21669,
148.0315, 30.34906264,
284.7158, 12.22117085,
220.1695, 4.284931111,
291.8024, 2.184704167
};
REAL8 sa[SDNUM];
/*
** delta long = lampl * COS (lphase - arg) in seconds of arc
** delta rad = rampl * COS (rphase - arg) in ninth place of log
** arg = j * sa (k) + i * ma (this planet)
** ma = mean anomaly
** sa = mean anomaly of disturbing planet, where this
** is taken from the aproximate value in sa[]
** For the COS (phase - arg) it is good enough to compute
** with 32 bit reals, because ampl and phase have only
** four to five significant digits.
** While saving constant space, it is costing execution time due
** to float/double conversions.
**
** In basic, all correction terms for sun, mercury, venus and mars
** were contained in one array K(0..142,0..6); Nk(N,0) contained
** the index of the first term of planet N and Nk(N,1) the number
** of terms for this planet. Here, we use a 0 in the first column
** kor.j to indicate the end of the table for a planet.
** K(*) was a basic INTEGER array, therefore the amplitudes and phases
** had to be expressed as
** K(i,2) = ampl. of longitude in 0.001 seconds of arc
** K(i,3) = phase of longitude in 0.01 degrees
** K(i,4) = ampl. of radius in 9th place of log
** K(i,5) = phase of radius in 0.01 degrees.
** Here we have converted the amplitude of long. to seconds of arc
** and the phases to degrees.
*/
struct kor ARR earthkor[86+1] = { /* 11-jul-88 all terms to 0.020" long */
/* j i lampl lphase rampl rphase k */
-1, 1, 0.013, 243, 28, 335, 8, /* mercury */
-1, 3, 0.015, 357, 18, 267, 8,
-1, 4, 0.023, 326, 5, 239, 8,
-1, 0, 0.075, 296.6, 94, 205.0, 0, /* venus */
-1, 1, 4.838, 299.10, 2359, 209.08, 0,
-1, 2, 0.074, 207.9, 69, 348.5, 0,
-1, 3, 0.009, 249, 16, 330, 0,
-2, 1, .116, 148.90, 160, 58.40, 0,
-2, 2, 5.526, 148.31, 6842, 58.32, 0,
-2, 3, 2.497, 315.94, 869, 226.70, 0,
-2, 4, 0.044, 311.4, 52, 38.8, 0,
-3, 2, 0.013, 176, 21, 90, 0,
-3, 3, .666, 177.71, 1045, 87.57, 0,
-3, 4, 1.559, -14.75, 1497, 255.25, 0,
-3, 5, 1.024, 318.15, 194, 49.50, 0,
-3, 6, 0.017, 315, 19, 43, 0,
-4, 4, .210, 206.20, 376, 116.28, 0,
-4, 5, .144, 195.40, 196, 105.20, 0,
-4, 6, .152, -16.20, 94, 254.80, 0,
-5, 5, 0.084, 235.6, 163, 145.4, 0,
-5, 6, 0.037, 221.8, 59, 132.2, 0,
-5, 7, .123, 195.30, 141, 105.40, 0,
-5, 8, .154, -.40, 26, 270.00, 0,
-6, 6, 0.038, 264.1, 80, 174.3, 0,
-6, 7, 0.014, 253, 25, 164, 0,
-6, 8, 0.01, 230, 14, 135, 0,
-6, 9, 0.014, 12, 12, 284, 0,
-7, 7, 0.020, 294, 42, 203.5, 0,
-7, 8, 0.006, 279, 12, 194, 0,
-8, 8, 0.011, 322, 24, 234, 0,
-8, 12, 0.042, 259.2, 44, 169.7, 0,
-8, 14, 0.032, 48.8, 33, 138.7, 0,
-9, 9, 0.006, 351, 13, 261, 0,
1, -1, .273, 217.70, 150, 127.70, 1, /* mars */
1, 0, 0.048, 260.3, 28, 347, 1,
2, -3, 0.041, 346, 52, 255.4, 1,
2, -2, 2.043, 343.89, 2057, 253.83, 1,
2, -1, 1.770, 200.40, 151, 295.00, 1,
2, 0, 0.028, 148, 31, 234.3, 1,
3, -3, .129, 294.20, 168, 203.50, 1,
3, -2, .425, -21.12, 215, 249.00, 1,
4, -4, 0.034, 71, 49, 339.7, 1,
4, -3, .500, 105.18, 478, 15.17, 1,
4, -2, .585, -25.94, 105, 65.90, 1,
5, -4, 0.085, 54.6, 107, 324.6, 1,
5, -3, .204, 100.80, 89, 11.00, 1,
6, -5, 0.020, 186, 30, 95.7, 1,
6, -4, .154, 227.40, 139, 137.30, 1,
6, -3, .101, 96.30, 27, 188.00, 1,
7, -5, 0.049, 176.5, 60, 86.2, 1,
7, -4, .106, 222.70, 38, 132.90, 1,
8, -5, 0.052, 348.9, 45, 259.7, 1,
8, -4, 0.021, 215.2, 8, 310, 1,
8, -6, 0.010, 307, 15, 217, 1,
9, -6, 0.028, 298, 34, 208.1, 1,
9, -5, 0.062, 346, 17, 257, 1,
10, -6, 0.019, 111, 15, 23, 1,
11, -7, 0.017, 59, 20, 330, 1,
11, -6, 0.044, 105.9, 9, 21, 1,
13, -8, 0.013, 184, 15, 94, 1,
13, -7, 0.045, 227.8, 5, 143, 1,
15, -9, 0.021, 309, 22, 220, 1,
17, -9, 0.026, 113, 0, 0, 1,
1, -2, .163, 198.60, 208, 112.00, 2, /* jupiter */
1, -1, 7.208, 179.53, 7067, 89.55, 2,
1, 0, 2.600, 263.22, 244, -21.40, 2,
1, 1, 0.073, 276.3, 80, 6.5, 2,
2, -3, 0.069, 80.8, 103, 350.5, 2,
2, -2, 2.731, 87.15, 4026, -2.89, 2,
2, -1, 1.610, 109.49, 1459, 19.47, 2,
2, 0, 0.073, 252.6, 8, 263, 2,
3, -3, .164, 170.50, 281, 81.20, 2,
3, -2, .556, 82.65, 803, -7.44, 2,
3, -1, .210, 98.50, 174, 8.60, 2,
4, -4, 0.016, 259, 29, 170, 2,
4, -3, 0.044, 168.2, 74, 79.9, 2,
4, -2, 0.080, 77.7, 113, 347.7, 2,
4, -1, 0.023, 93, 17, 3, 2,
5, -2, 0.009, 71, 14, 343, 2,
1, -2, 0.011, 105, 15, 11, 3, /* saturn */
1, -1, .419, 100.58, 429, 10.60, 3,
1, 0, .320, 269.46, 8, -7.00, 3,
2, -2, .108, 290.60, 162, 200.60, 3,
2, -1, .112, 293.60, 112, 203.10, 3,
3, -2, 0.021, 289, 32, 200.1, 3,
3, -1, 0.017, 291, 17, 201, 3,
ENDMARK
};
struct kor ARR mercurykor[24+1] = {
1, -1, .711, 35.47, 491, 305.28, 4,
2, -3, .552, 161.15, 712, 71.12, 4,
2, -2, 2.100, 161.15, 2370, 71.19, 4,
2, -1, 3.724, 160.69, 899, 70.49, 4,
2, 0, .729, 159.76, 763, 250.00, 4,
3, -3, .431, 105.37, 541, 15.53, 4,
3, -2, 1.329, 104.78, 1157, 14.84, 4,
3, -1, .539, 278.95, 14, 282.00, 4,
4, -2, .484, 226.40, 234, 136.02, 4,
5, -4, .685, -10.43, 849, 259.51, 4,
5, -3, 2.810, -10.14, 2954, 259.92, 4,
5, -2, 7.356, -12.22, 282, 255.43, 4,
5, -1, 1.471, -12.30, 1550, 77.75, 4,
5, 0, .375, -12.29, 472, 77.70, 4,
2, -1, .443, 218.48, 256, 128.36, 5,
4, -2, .374, 151.81, 397, 61.63, 5,
4, -1, .808, 145.93, 13, 35.00, 5,
1, -1, .697, 181.07, 708, 91.38, 6,
1, 0, .574, 236.72, 75, 265.40, 6,
2, -2, .938, 36.98, 1185, 306.97, 6,
2, -1, 3.275, 37.00, 3268, 306.99, 6,
2, 0, .499, 31.91, 371, 126.90, 6,
3, -1, .353, 25.84, 347, 295.76, 6,
2, -1, .380, 239.87, 0, 0, 7,
ENDMARK
};
struct kor ARR venuskor[22+1] = {
-1, 2, .264, -19.20, 175, 251.10, 8,
-2, 5, .361, 167.68, 55, 77.20, 8,
1, -1, 4.889, 119.11, 2246, 29.11, 9,
2, -2, 11.261, 148.23, 9772, 58.21, 9,
3, -3, 7.128, -2.57, 8271, 267.42, 9,
3, -2, 3.446, 135.91, 737, 47.37, 9,
4, -4, 1.034, 26.54, 1426, 296.49, 9,
4, -3, .677, 165.32, 445, 75.70, 9,
5, -5, .330, 56.88, 510, -33.36, 9,
5, -4, 1.575, 193.93, 1572, 104.21, 9,
5, -3, 1.439, 138.08, 162, 229.90, 9,
6, -6, .143, 84.40, 236, -5.80, 9,
6, -5, .205, 44.20, 256, 314.20, 9,
6, -4, .176, 164.30, 70, 75.70, 9,
8, -5, .231, 180.00, 25, 75.00, 9,
3, -2, .673, 221.62, 717, 131.60, 10,
3, -1, 1.208, 237.57, 29, 149.00, 10,
1, -1, 2.966, 208.09, 2991, 118.09, 11,
1, 0, 1.563, 268.31, 91, -7.60, 11,
2, -2, .889, 145.16, 1335, 55.17, 11,
2, -1, .480, 171.01, 464, 80.95, 11,
3, -2, .169, 144.20, 250, 54.00, 11,
ENDMARK
};
struct kor ARR marskor[62+1] = {
-1, 1, .115, 65.84, 684, 156.14, 12,
-1, 2, .623, 246.03, 812, 155.77, 12,
-1, 3, 6.368, 57.60, 556, -32.06, 12,
-1, 4, .588, 57.24, 616, 147.28, 12,
-2, 5, .138, 39.18, 157, 309.39, 12,
-2, 6, .459, 217.58, 82, 128.10, 12,
-1, -1, .106, 33.60, 141, 303.45, 13,
-1, 0, .873, 34.34, 1112, 304.05, 13,
-1, 1, 8.559, 35.10, 6947, 304.45, 13,
-1, 2, 13.966, 20.50, 2875, 113.20, 13,
-1, 3, 1.487, 22.18, 1619, 112.38, 13,
-1, 4, .175, 22.46, 225, 112.15, 13,
-2, 2, .150, 18.96, 484, 266.42, 13,
-2, 3, 7.355, 158.64, 6412, 68.62, 13,
-2, 4, 4.905, 154.09, 1985, 244.70, 13,
-2, 5, .489, 154.39, 543, 244.50, 13,
-3, 3, .216, 111.06, 389, 21.10, 13,
-3, 4, .355, 110.64, 587, 19.17, 13,
-3, 5, 2.641, 280.58, 2038, 190.60, 13,
-3, 6, .970, 276.06, 587, 6.75, 13,
-3, 7, .100, 276.20, 116, 6.40, 13,
-4, 5, .152, 232.48, 259, 142.60, 13,
-4, 6, .264, 230.47, 387, 139.75, 13,
-4, 7, 1.156, 41.64, 749, 312.67, 13,
-4, 8, .259, 37.92, 205, 128.80, 13,
-5, 8, .172, -8.99, 234, 260.70, 13,
-5, 9, .575, 164.48, 308, 74.60, 13,
-6, 10, .115, 113.70, 145, 23.53, 13,
-6, 11, .363, 285.69, 144, 196.00, 13,
-7, 13, .353, 48.83, 85, 319.10, 13,
-8, 15, 1.553, 170.14, 110, 81.00, 13,
-8, 16, .148, 170.74, 154, 259.94, 13,
-9, 17, .193, 293.70, 23, 22.80, 13,
1, -3, .382, 46.48, 521, 316.25, 14,
1, -2, 3.144, 46.78, 3894, 316.39, 14,
1, -1, 25.384, 48.96, 23116, 318.87, 14,
1, 0, 3.732, -17.62, 1525, 117.81, 14,
1, 1, .474, -34.60, 531, 59.67, 14,
2, -4, .265, 192.88, 396, 103.12, 14,
2, -3, 2.108, 192.72, 3042, 102.89, 14,
2, -2, 16.035, 191.90, 22144, 101.99, 14,
2, -1, 21.869, 188.35, 16624, 98.33, 14,
2, 0, 1.461, 189.66, 1478, 279.04, 14,
2, 1, .167, 191.04, 224, 280.81, 14,
3, -4, .206, 167.11, 338, 76.13, 14,
3, -3, 1.309, 168.27, 2141, 76.24, 14,
3, -2, 2.607, 228.41, 3437, 139.74, 14,
3, -1, 3.174, 207.20, 1915, 115.83, 14,
3, 0, .232, 207.78, 240, 298.06, 14,
4, -4, .178, 127.25, 322, 36.16, 14,
4, -3, .241, 200.69, 389, 110.02, 14,
4, -2, .330, 267.57, 413, 179.86, 14,
4, -1, .416, 221.88, 184, 128.17, 14,
1, -2, .155, -38.20, 191, 231.58, 15,
1, -1, 1.351, -34.10, 1345, 235.85, 15,
1, 0, .884, 288.05, 111, 39.90, 15,
1, 1, .132, 284.88, 144, 15.67, 15,
2, -2, .620, 35.15, 869, 305.30, 15,
2, -1, 1.768, 32.50, 1661, 302.51, 15,
2, 0, .125, 18.73, 103, 119.90, 15,
3, -2, .141, 47.59, 199, 318.06, 15,
3, -1, .281, 40.95, 248, 310.75, 15,
ENDMARK
};
struct m45dat m45[NUM_MOON_CORR] = {
/* l, l', F, D, Long, Lat, Par),*/
{ 0, 0, 0, 4, 13.902, 14.06, 0.2607},
{ 0, 0, 0, 2, 2369.912, 2373.36, 28.2333},
{ 1, 0, 0, 4, 1.979, 6.98, 0.0433},
{ 1, 0, 0, 2, 191.953, 192.72, 3.0861},
{ 1, 0, 0, 0, 22639.500, 22609.1, 186.5398},
{ 1, 0, 0, -2, -4586.465, -4578.13, 34.3117},
{ 1, 0, 0, -4, -38.428, -38.64, 0.6008},
{ 1, 0, 0, -6, -0.393, -1.43, 0.0086},
{ 0, 1, 0, 4, -0.289, -1.59, -0.0053},
{ 0, 1, 0, 2, -24.420, -25.10, -0.3000},
{ 0, 1, 0, 0, -668.146, -126.98, -0.3997},
{ 0, 1, 0, -2, -165.145, -165.06, 1.9178},
{ 0, 1, 0, -4, -1.877, -6.46, 0.0339},
{ 0, 0, 0, 3, 0.403, -4.01, 0.0023},
{ 0, 0, 0, 1, -125.154, -112.79, -0.9781},
{ 2, 0, 0, 4, 0.213, 1.02, 0.0054},
{ 2, 0, 0, 2, 14.387, 14.78, 0.2833},
{ 2, 0, 0, 0, 769.016, 767.96, 10.1657},
{ 2, 0, 0, -2, -211.656, -152.53, -0.3039},
{ 2, 0, 0, -4, -30.773, -34.07, 0.3722},
{ 2, 0, 0, -6, -0.570, -1.40, 0.0109},
{ 1, 1, 0, 2, -2.921, -11.75, -0.0484},
{ 1, 1, 0, 0, -109.673, -115.18, -0.9490},
{ 1, 1, 0, -2, -205.962, -182.36, 1.4437},
{ 1, 1, 0, -4, -4.391, -9.66, 0.0673},
{ 1, -1, 0, 4, 0.283, 1.53, 0.0060},
{ 1, -1, 0, 2, 14.577, 31.70, 0.2302},
{ 1, -1, 0, 0, 147.687, 138.76, 1.1528},
{ 1, -1, 0, -2, 28.475, 23.59, -0.2257},
{ 1, -1, 0, -4, 0.636, 2.27, -0.0102},
{ 0, 2, 0, 2, -0.189, -1.68, -0.0028},
{ 0, 2, 0, 0, -7.486, -0.66, -0.0086},
{ 0, 2, 0, -2, -8.096, -16.35, 0.0918},
{ 0, 0, 2, 2, -5.741, -0.04, -0.0009},
{ 0, 0, 2, 0, -411.608, -0.2, -0.0124},
{ 0, 0, 2, -2, -55.173, -52.14, -0.1052},
{ 0, 0, 2, -4, 0.025, -1.67, 0.0031},
{ 1, 0, 0, 1, -8.466, -13.51, -0.1093},
{ 1, 0, 0, -1, 18.609, 3.59, 0.0118},
{ 1, 0, 0, -3, 3.215, 5.44, -0.0386},
{ 0, 1, 0, 1, 18.023, 17.93, 0.1494},
{ 0, 1, 0, -1, 0.560, 0.32, -0.0037},
{ 3, 0, 0, 2, 1.060, 2.96, 0.0243},
{ 3, 0, 0, 0, 36.124, 50.64, 0.6215},
{ 3, 0, 0, -2, -13.193, -16.40, -0.1187},
{ 3, 0, 0, -4, -1.187, -0.74, 0.0074},
{ 3, 0, 0, -6, -0.293, -0.31, 0.0046},
{ 2, 1, 0, 2, -0.290, -1.45, -0.0051},
{ 2, 1, 0, 0, -7.649, -10.56, -0.1038},
{ 2, 1, 0, -2, -8.627, -7.59, -0.0192},
{ 2, 1, 0, -4, -2.740, -2.54, 0.0324},
{ 2, -1, 0, 2, 1.181, 3.32, 0.0213},
{ 2, -1, 0, 0, 9.703, 11.67, 0.1268},
{ 2, -1, 0, -2, -2.494, -1.17, -0.0017},
{ 2, -1, 0, -4, 0.360, 0.20, -0.0043},
{ 1, 2, 0, 0, -1.167, -1.25, -0.0106},
{ 1, 2, 0, -2, -7.412, -6.12, 0.0484},
{ 1, 2, 0, -4, -0.311, -0.65, 0.0044},
{ 1, -2, 0, 2, 0.757, 1.82, 0.0112},
{ 1, -2, 0, 0, 2.580, 2.32, 0.0196},
{ 1, -2, 0, -2, 2.533, 2.40, -0.0212},
{ 0, 3, 0, -2, -0.344, -0.57, 0.0036},
{ 1, 0, 2, 2, -0.992, -0.02, 0},
{ 1, 0, 2, 0, -45.099, -0.02, -0.0010},
{ 1, 0, 2, -2, -0.179, -9.52, -0.0833},
{ 1, 0, -2, 2, -6.382, -3.37, -0.0481},
{ 1, 0, -2, 0, 39.528, 85.13, -0.7136},
{ 1, 0, -2, -2, 9.366, 0.71, -0.0112},
{ 0, 1, 2, 0, 0.415, 0.10, 0.0013},
{ 0, 1, 2, -2, -2.152, -2.26, -0.0066},
{ 0, 1, -2, 2, -1.440, -1.30, 0.0014},
{ 0, 1, -2, -2, 0.384, 0.0, 0.0},
{ 2, 0, 0, 1, -0.586, -1.20, -0.0100},
{ 2, 0, 0, -1, 1.750, 2.01, 0.0155},
{ 2, 0, 0, -3, 1.225, 0.91, -0.0088},
{ 1, 1, 0, 1, 1.267, 1.52, 0.0164},
{ 1, -1, 0, -1, -1.089, 0.55, 0},
{ 0, 0, 2, -1, 0.584, 8.84, 0.0071},
{ 4, 0, 0, 0, 1.938, 3.60, 0.0401},
{ 4, 0, 0, -2, -0.952, -1.58, -0.0130},
{ 3, 1, 0, 0, -0.551, 0.94, -0.0097},
{ 3, 1, 0, -2, -0.482, -0.57, -0.0045},
{ 3, -1, 0, 0, 0.681, 0.96, 0.0115},
{ 2, 0, 2, 0, -3.996, 0, 0.0004},
{ 2, 0, 2, -2, 0.557, -0.75, -0.0090},
{ 2, 0, -2, 2, -0.459, -0.38, -0.0053},
{ 2, 0, -2, 0, -1.298, 0.74, 0.0004},
{ 2, 0, -2, -2, 0.538, 1.14, -0.0141},
{ 1, 1, -2, -2, 0.426, 0.07, -0.0006},
{ 1, -1, 2, 0, -0.304, 0.03, 0.0003},
{ 1, -1, -2, 2, -0.372, -0.19, -0.0027},
{ 0, 0, 4, 0, 0.418, 0, 0},
{ 2, -1, 0, -1, -0.352, -0.37, -0.0028}
};
# if MOON_TEST_CORR
/* moon additional correction terms */
struct m5dat {
REAL8 lng;
int i0,i1,i2,i3;
} m5[] = {
/* lng, l, l', F, D, */
0.127, 0, 0, 0, 6,
-0.151, 0, 2, 0, -4,
-0.085, 0, 0, 2, 4,
0.150, 0, 1, 0, 3,
-0.091, 2, 1, 0, -6,
-0.103, 0, 3, 0, 0,
-0.301, 1, 0, 2, -4,
0.202, 1, 0, -2, -4,
0.137, 1, 1, 0, -1,
0.233, 1, 1, 0, -3,
-0.122, 1, -1, 0, 1,
-0.276, 1, -1, 0, -3,
0.255, 0, 0, 2, 1,
0.254, 0, 0, 2, -3,
-0.100, 3, 1, 0, -4,
-0.183, 3, -1, 0, -2,
-0.297, 2, 2, 0, -2,
-0.161, 2, 2, 0, -4,
0.197, 2, -2, 0, 0,
0.254, 2, -2, 0, -2,
-0.250, 1, 3, 0, -2,
-0.123, 2, 0, 2, 2,
0.173, 2, 0, -2, -4,
0.263, 1, 1, 2, 0,
0.130, 3, 0, 0, -1,
0.113, 5, 0, 0, 0,
0.092, 3, 0, 2, -2,
0, 99, 0, 0, 0 /* end mark */
};
# endif /* MOON_TEST_CORR */
/* solution of the Kepler equation, return rad */
/* t = mean anomaly in degrees */
/* ex = excentricity of orbit */
/* err = maximum error in degrees */
REAL8 fnu(t, ex, err)
REAL8 t;
REAL8 ex;
REAL8 err;
{
REAL8 u0, delta;
t *= DEGTORAD;
u0 = t;
err *= DEGTORAD;
delta = 1;
while (ABS8(delta) >= err) {
delta = (t + ex * SIN8(u0) - u0) / (1 - ex * COS8(u0));
u0 += delta;
}
return u0;
}
/* x MOD 360.0, so that x in 0..360 */
REAL8 smod8360(x)
REAL8 x;
{
while (x >= 360.0)
x -= 360.0;
while (x < 0.0)
x += 360.0;
return x;
}
/* x MOD 360.0, so that x in 0..360 */
REAL8 mod8360(x)
REAL8 x;
{
if (x >= 0 && x < 360.0)
return x;
return x - 360.0 * RFloor(x / 360.0);
}
/* a - b on a 360 degree circle, result -180..180 */
REAL8 diff8360(a, b)
REAL8 a;
REAL8 b;
{
REAL8 d;
d = a - b;
if (d >= 180.0)
return d - 360.0;
if (d < -180.0)
return d + 360.0;
return d;
}
REAL8 test_near_zero(x)
REAL8 x;
{
if (ABS8(x) >= NEAR_ZERO)
return x;
if (x < 0)
return -NEAR_ZERO;
return NEAR_ZERO;
}
/*
** p points to memory filled with long values; for
** each of the values the seqeuence of the four bytes
** has to be reversed, to translate HP-UX and VAX
** ordering to MSDOS/Turboc ordering
*/
void longreorder(p, n)
UCHAR *p;
int n;
{
int i;
unsigned char c0, c1, c2, c3;
static int orderinit = 0;
unsigned short test;
if (!orderinit) {
test = 0x0001;
orderinit = (*(unsigned char *)(&test)) ? 1 : -1;
}
if (orderinit < 0)
return;
for (i = 0; i < n; i += 4, p += 4) {
c0 = *p;
c1 = *(p + 1);
c2 = *(p + 2);
c3 = *(p + 3);
*p = c3;
*(p + 1) = c2;
*(p + 2) = c1;
*(p + 3) = c0;
}
}
/*****************************************************
$Header: deltat.c,v 1.10 93/01/27 14:37:06 alois Exp $
deltat.c
deltat(t): returns delta t (in julian days) from universal time t
is included by users
ET = UT + deltat
---------------------------------------------------------------
| Copyright Astrodienst Zurich AG and Alois Treindl, 1989. |
| The use of this source code is subject to regulations made |
| by Astrodienst Zurich. The code is NOT in the public domain.|
| |
| This copyright notice must not be changed or removed |
| by any user of this program. |
---------------------------------------------------------------
******************************************************/
double deltat(jd_ad)
double jd_ad;
{
static short int ARR dt[] = { /* in centiseconds */
/*
** dt from 1637 to 2000, as tabulated in A.E.
** the values 1620 - 1636 are not taken, as they fit
** badly the parabola 25.5 t*t for the next range. The
** best crossing point to switch to the parabola is
** 1637, where we have fitted the value for continuity
*/
6780, 6500, 6300,
6200, 6000, 5800, 5700, 5500,
5400, 5300, 5100, 5000, 4900,
4800, 4700, 4600, 4500, 4400,
4300, 4200, 4100, 4000, 3800, /* 1655 - 59 */
3700, 3600, 3500, 3400, 3300,
3200, 3100, 3000, 2800, 2700,
2600, 2500, 2400, 2300, 2200,
2100, 2000, 1900, 1800, 1700,
1600, 1500, 1400, 1400, 1300,
1200, 1200, 1100, 1100, 1000,
1000, 1000, 900, 900, 900,
900, 900, 900, 900, 900,
900, 900, 900, 900, 900, /* 1700 - 1704 */
900, 900, 900, 1000, 1000,
1000, 1000, 1000, 1000, 1000,
1000, 1000, 1100, 1100, 1100,
1100, 1100, 1100, 1100, 1100,
1100, 1100, 1100, 1100, 1100,
1100, 1100, 1100, 1100, 1200, /* 1730 - 1734 */
1200, 1200, 1200, 1200, 1200,
1200, 1200, 1200, 1200, 1300,
1300, 1300, 1300, 1300, 1300,
1300, 1400, 1400, 1400, 1400,
1400, 1400, 1400, 1500, 1500,
1500, 1500, 1500, 1500, 1500, /* 1760 - 1764 */
1600, 1600, 1600, 1600, 1600,
1600, 1600, 1600, 1600, 1600,
1700, 1700, 1700, 1700, 1700,
1700, 1700, 1700, 1700, 1700,
1700, 1700, 1700, 1700, 1700,
1700, 1700, 1600, 1600, 1600, /* 1790 - 1794 */
1600, 1500, 1500, 1400, 1400,
1370, 1340, 1310, 1290, 1270, /* 1800 - 1804 */
1260, 1250, 1250, 1250, 1250,
1250, 1250, 1250, 1250, 1250,
1250, 1250, 1240, 1230, 1220,
1200, 1170, 1140, 1110, 1060,
1020, 960, 910, 860, 800,
750, 700, 660, 630, 600, /* 1830 - 1834 */
580, 570, 560, 560, 560,
570, 580, 590, 610, 620,
630, 650, 660, 680, 690,
710, 720, 730, 740, 750,
760, 770, 770, 780, 780,
788, 782, 754, 697, 640, /* 1860 - 1864 */
602, 541, 410, 292, 182,
161, 10, -102, -128, -269,
-324, -364, -454, -471, -511,
-540, -542, -520, -546, -546,
-579, -563, -564, -580, -566,
-587, -601, -619, -664, -644, /* 1890 - 1894 */
-647, -609, -576, -466, -374,
-272, -154, -2, 124, 264,
386, 537, 614, 775, 913,
1046, 1153, 1336, 1465, 1601,
1720, 1824, 1906, 2025, 2095,
2116, 2225, 2241, 2303, 2349, /* 1920 - 1924 */
2362, 2386, 2449, 2434, 2408,
2402, 2400, 2387, 2395, 2386,
2393, 2373, 2392, 2396, 2402,
2433, 2483, 2530, 2570, 2624,
2677, 2728, 2778, 2825, 2871,
2915, 2957, 2997, 3036, 3072, /* 1950 - 1954 */
3107, 3135, 3168, 3218, 3268,
3315, 3359, 3400, 3447, 3503,
3573, 3654, 3743, 3829, 3920,
4018, 4117, 4223, 4337, 4449,
4548, 4646, 4752, 4853, 4959,
5054, 5138, 5217, 5296, 5379, /* 1980 - 1984 */
5434, 5487, 5532, 5582, 5630, /* 1985 - 89 from AE 1991 */
5686, 5757, 5900, 5900, 6000, /* AE 1993 and extrapol */
6050, 6100, 6150, 6200, 6250, /* 1995 - 1999 */
6300}; /* 2000 */
double yr, cy, delta;
long iyr, i;
yr = (jd_ad + 18262) / 365.25 + 100.0; /* year relative 1800 */
cy = yr / 100;
iyr = (long) (RFloor(yr) + 1800); /* truncated to integer, rel 0 */
#if TIDAL_26 /* Stephenson formula only when 26" tidal
term in lunar motion */
if (iyr >= 1637 && iyr < 2000) {
i = iyr - 1637;
delta = dt[i] * 0.01 + (dt[i+1] - dt[i]) * (yr - RFloor(yr)) * 0.01;
} else if (iyr >= 2000) { /* parabola, fitted at value[2000] */
delta = 25.5 * cy * cy - 25.5 * 4 + 63.00;
} else if (iyr >= 948) { /* from 948 - 1637 use parabola */
delta = 25.5 * cy * cy;
} else { /* before 984 use other parabola */
delta = 1361.7 + 320 * cy + 44.3 * cy * cy; /* fits at 948 */
}
#else /* use Clemence value + 5 sec before 1690, new dt afterwards */
cy -= 1; /* epoch 1900 */
if (iyr >= 1690 && iyr < 2000) {
i = iyr - 1637;
delta = dt[i] * 0.01 + (dt[i+1] - dt[i]) * (yr - RFloor(yr)) * 0.01;
} else if (iyr >= 2000) { /* parabola, fitted at value[2000] */
delta = 29.949 * cy * cy - 29.949 * 4 + 63.0;
} else {
delta = 5 + 24.349 + 72.3165 * cy + 29.949 * cy * cy; /* fits at 1690 */
}
#endif
return delta / 86400.0;
}
/*******************************************
$Header: d2l.c,v 1.9 91/11/16 16:24:20 alois Exp $
double to long with rounding, no overflow check
*************************************/
long d2l(x)
double x;
{
if (x >=0)
return ((long) (x + 0.5));
else
return (-(long) (0.5 - x));
}
/*
* $Header$
*
* A collection of useful functions for centisec calculations.
---------------------------------------------------------------
| Copyright Astrodienst Zurich AG and Alois Treindl, 1991. |
| The use of this source code is subject to regulations made |
| by Astrodienst Zurich. The code is NOT in the public domain.|
| |
| This copyright notice must not be changed or removed |
| by any user of this program. |
---------------------------------------------------------------
*******************************************************/
double degnorm(p)
double p;
{
if (p < 0)
do {
p += 360.0;
} while (p < 0);
else if (p >= 360.0)
do {
p -= 360.0;
} while (p >= 360.0);
return (p);
}
/*********************************************************
$Header: julday.c,v 1.9 91/11/16 16:25:06 alois Exp $
*********************************************************/
/*
** This function returns the absolute Julian day number (JD)
** for a given calendar date.
** The aruguments are a calendar date: day, month, year as integers,
** hour as double with decimal fraction.
** If gregflag = 1, Gregorian calendar is assumed, gregflag = 0
** Julian calendar is assumed.
**
** The Julian day number is system of numbering all days continously
** within the time range of known human history. It should be familiar
** for every astrological or astronomical programmer. The time variable
** in astronomical theories is usually expressed in Julian days or
** Julian centuries (36525 days per century) relative to some start day;
** the start day is called 'the epoch'.
** The Julian day number is a double representing the number of
** days since JD = 0.0 on 1 Jan -4712, 12:00 noon.
** Midnight has always a JD with fraction .5, because traditionally
** the astronomical day started at noon.
**
** NOTE: The Julian day number is named after the monk Julianus. It must
** not be confused with the Julian calendar system, which is named after
** Julius Cesar, the Roman politician who introduced this calendar.
**
** Original author: Marc Pottenger, Los Angeles.
** with bug fix for year < -4711 15-aug-88 by Alois Treindl
**
** References: Oliver Montenbruck, Grundlagen der Ephemeridenrechnung,
** Verlag Sterne und Weltraum (1987), p.49 ff
**
** related functions: revjul() reverse Julian day number: compute the
** calendar date from a given JD
*/
double julday(month, day, year, hour, gregflag)
int month;
int day;
int year;
double hour;
int gregflag;
{
double jd, u, u0, u1, u2;
u = year;
if (month < 3)
u -=1;
u0 = u + 4712.0;
u1 = month + 1.0;
if (u1 < 4)
u1 += 12.0;
jd = RFloor(u0*365.25)
+ RFloor(30.6*u1+0.000001)
+ day + hour/24.0 - 63.5;
if (gregflag) {
u2 = RFloor(ABS8(u) / 100) - RFloor(ABS8(u) / 400);
if (u < 0.0)
u2 = -u2;
jd = jd - u2 + 2;
if ((u < 0.0) && (u/100 == RFloor(u/100)) && (u/400 != RFloor(u/400)))
jd -= 1;
}
return jd;
}
/*********************************************************
$Header: revjul.c,v 1.9 91/11/16 16:25:37 alois Exp $
*********************************************************/
/*
** revjul() is the inverse function to julday(), see the description there.
** Arguments are julian day number, calendar flag (0=julian, 1=gregorian)
** return values are the calendar day, month, year and the hour of
** the day with decimal fraction (0 .. 23.999999).
**
** Original author Mark Pottenger, Los Angeles.
** with bug fix for year < -4711 16-aug-88 Alois Treindl
*/
void revjul(jd, gregflag, jmon, jday, jyear, jut)
double jd;
int gregflag;
int *jmon;
int *jday;
int *jyear;
double *jut;
{
double u0, u1, u2, u3, u4;
u0 = jd + 32082.5;
if (gregflag) {
u1 = u0 + RFloor(u0/36525.0) - RFloor(u0/146100.0) - 38.0;
if (jd >= 1830691.5) u1 +=1;
u0 = u0 + RFloor(u1/36525.0) - RFloor(u1/146100.0) - 38.0;
}
u2 = RFloor(u0 + 123.0);
u3 = RFloor((u2 - 122.2) / 365.25);
u4 = RFloor((u2 - RFloor(365.25 * u3)) / 30.6001);
*jmon = (int)(u4-1.0);
if (*jmon > 12)
*jmon -= 12;
*jday = (int)(u2 - RFloor(365.25 * u3) - RFloor(30.6001 * u4));
*jyear = (int)(u3 + RFloor((u4 - 1.9999) / 12.0) - 4800.0);
*jut = (jd - RFloor(jd + 0.5) + 0.5) * 24.0;
}
#endif /* PLACALC */
/* placalc2.c */